The impact of ZIKV infection on gene expression in neural cells over time

Zika virus (ZIKV) outbreak caused one of the most significant medical emergencies in the Americas due to associated microcephaly in newborns. To evaluate the impact of ZIKV infection on neuronal cells over time, we retrieved gene expression data from several ZIKV-infected samples obtained at different time point post-infection (pi). Differential gene expression analysis was applied at each time point, with more differentially expressed genes (DEG) identified at 72h pi. There were 5 DEGs (PLA2G2F, TMEM71, PKD1L2, UBD, and TNFAIP3 genes) across all timepoints, which clearly distinguished between infected and healthy samples. The highest expression levels of all five genes were identified at 72h pi. Taken together, our results indicate that ZIKV infection greatly impacts human neural cells at early times of infection, with peak perturbation observed at 72h pi. Our analysis revealed that all five DEGs, in samples of ZIKV-infected human neural stem cells, remained highly upregulated across the timepoints evaluated. Moreover, despite the pronounced inflammatory host response observed throughout infection, the impact of ZIKV is variable over time. Finally, the five DEGs identified herein play prominent roles in infection, and could serve to guide future investigations into virus-host interaction, as well as constitute targets for therapeutic drug development.


Introduction
The Zika Virus (ZIKV) outbreak in the Americas captured the world's attention in 2016, leading to one of the most important medical emergencies in the last decade [1].ZIKV is a mosquito-borne Flavivirus belonging to the Flaviviridae family and is closely related to other viruses (e.g., West Nile Virus, Dengue virus, and Yellow fever virus) that cause arthropodborne diseases in humans [2].Generally, these diseases are transmitted to humans by females of infected mosquitoes belonging to the Aedes genus during blood-feeding [3].
In Brazil, the first confirmed case of ZIKV was reported in 2015 [4], followed by the swift spread of the virus [5].In 2016, more than 26 countries across the Americas and Europe reported cases of the disease [6].Most often (in over 80% of ZIKV infections in healthy individuals), disease presentation is an asymptomatic infection or characterized by mild symptoms [7][8][9].However, less than one year after the first reported case in Salvador, the capital of the state of Bahia-Brazil, increasing reports of Guillain-Barre ´syndrome (GBS) was noted [10].Furthermore, in the state of Pernambuco, prenatal ZIKV exposure was associated with an almost 20-fold increase in the incidence of microcephaly [11].Currently, ZIKV infection during pregnancy could lead to what is known as Congenital Zika Syndrome (CZS).CZS has been linked to brain abnormalities and/or microcephaly or neural tube defects, eye abnormalities, or consequent central nervous system dysfunction among fetuses or infants in the absence of other evidence of brain abnormalities or microcephaly [12].
Since the 2016 outbreak, the effects of ZIKV infection on gene expression in the host's brain cells have been extensively studied to elucidate the molecular mechanisms related to CZS neuropathogenesis.A previous study explored the effects of ZIKV infection on mRNA and microRNA (miRNA) expression in astrocytes, aiming to identify important miRNAs and pathways related to viral pathology [13].The effects of ZIKV infection on human neural stem cells (hNSCs) have also recently been studied [14].Further investigations established the importance of viral infection in hNSCs with respect to the neuropathogenesis of CZS [15][16][17], thereby corroborating the previous study's findings, by identifying several miRNAs related to neurodevelopment and oxidative stress.Other studies have also provided insights about ZIKV's effects on miRNA expression in hNSCs [18], highlighting miRNAs that mediate the suppression of gene networks related to cell cycle progression and stem cell maintenance.Nonetheless, the pathogenesis of ZIKV in neuronal cells has not been completely elucidated, and an mRNA-focused biosignature is still absent from the current literature.
In order to provide further insights into the effects of ZIKV infection on neuronal stem cells and potentially enhance the body of knowledge on CZS development, an integrated gene expression analysis was performed in publicly available gene expression datasets.A similar approach has been widely employed to identify important genes and molecular pathways with the aim of discovering potential biomarkers in several pathologies, such as endometrial cancer [19], Sickle cell disease [20], and HTLV-1 [21].Related studies have analyzed Aedes aegypti expression datasets to identify important genes and pathways related to mosquito infection by Dengue, Yellow fever, West Nile, and Zika viruses [22,23].The present investigation analyzed neuronal stem cell gene expression data among 39 samples of which 30 were ZIKV-infected samples collected at different times post-infection and 9 non infected samples.

Data acquisition
The flow chart for the study selection for this integrated analysis is shown in Fig 1 .We have performed a dataset search on the Gene Expression Omnibus (GEO) repository to assess data about the gene expression of neuronal stem cells infected by ZIKV [24].86 GEO datasets were founded; however, 77 records were excluded because the samples were not induced pluripotent stem cell-derived human neural stem cells, resulting in 9 datasets.Of these, seven were excluded because they do not have more than one timepoint or because they only have information on non-coding RNAs.Thus, only two datasets were included (GSE97872 and GSE157530).
The data from both datasets were downloaded from the GEO repository [24] using the GEOquery package [25].Given that these are public datasets, Ethics Committees authorizations are not required for this paper.A total of 39 samples from human neural stem cells were retrieved, comprising 9 uninfected samples (healthy controls), 6 ZIKV-infected samples collected at 24 hours post-infection (pi), 6 samples at 48h pi, 6 samples at 72h pi, 3 samples at 5 weeks pi, 3 samples at 8 weeks pi,3 samples at 11 weeks pi, and 3 samples at 14 weeks pi.The raw data were log2-transformed and batch-effect corrected using the sva package [26] and then subjected to downstream analysis.

Differentially expressed genes
Differential expression analysis was employed to compare gene expression profiles in ZIKVinfected hNSCs over time with healthy controls in order to identify differentially expressed genes (DEGs).This analysis was performed using lmFit function present on limma package [27] in R 4.1.0.The function fit a linear model assuming Gaussian probability distribution on data than, DEGs were identified using a log fold-change threshold of ±1, with a false discovery rate corrected p-value < 0.05.Comparisons of DEGs at different sample collection times were performed by constructing Venn diagrams using the VennDiagram package [28].The com-pareCluster package [29] was then used to scan the REACTOME pathway database [30] for the obtained DEGs to identify enriched pathways at each selected timepoint.All data and values are available in S1 Table.

Evaluation of sample heterogeneity over time
Variation among samples within and between timepoints was evaluated using the molecular degree of perturbation (MDP) package [31], applied to gene expression values following batch correction.All comparisons between timepoints were performed with the Kruskal-Wallis test, followed by Tukey's posttest.

Data visualization and dimensionality reduction
To evaluate sample clustering and classification over time, we performed one-sided, unsupervised Ward's method hierarchical cluster analysis [32], heatmaps [33] and principal component analysis (PCA) using the obtained gene expression values.This approach allowed for the visualization of sample dispersion across timepoints.
We then investigated the pathways associated with specific DEGs identified at each timepoint.This analysis revealed an intense inflammatory response via the enriched pathways Interferon alpha/beta signaling and Interferon Signaling at 48h pi, which persisted until 14 .At 24h pi, the enriched pathways Extracellular matrix organization, Integrin cell surface interactions, Immunoregulatory interactions between a Lymphoid and a non-Lymphoid cell, and Degradation of the extracellular matrix indicated a cell organization and interaction response.At 72h pi, an intensified immune response was evidenced via the enrichment of the pathways Interleukin-4 and Interleukin-13 signaling, NGF-stimulated transcription, Nuclear Events (kinase and transcription factor activation) and Antigen Presentation: Folding, assembly and peptide loading of class I MHC.At 5 weeks pi, a specific pattern of alterations in synapses and cell-signaling was identified via the enriched pathways Regulation of FZD by ubiquitination, Transmission across Chemical Synapses, RUNX2 regulates osteoblast differentiation Trafficking of AMPA receptors, Glutamate binding, activation of AMPA receptors and synaptic plasticity and RUNX2 regulates bone development.However, at 8 weeks pi, an unspecific pattern related to the enrichment of pathways ATF4 activates genes in response to endoplasmic reticulum stress, PERK regulates gene expression, Unfolded Protein Response (UPR) and Response of EIF2AK4 (GCN2) to amino acid deficiency was identified.At 11 weeks pi, just one pathway was enriched: Ovarian tumor domain proteases (Fig 3).
We assessed the sample variations using degree of molecular perturbation and Principal Component Analysis (PCA) in the ZIKV samples to investigate variation across timepoints.We observed variation in ZIKV-infection-associated gene expression at all pi timepoints investigated.The highest variation was identified at 72h pi, significantly greater than at all other timepoints except week 5 (Fig 4A and 4B), suggesting that ZIKV may exert a stronger impact 72 hours after being bitten by an infected mosquito, with a progressive reduction in intensity over time.
The five genes that were significantly differentially expressed at all timepoints were investigated.All were found to be upregulated across the timepoints, with the highest log2 foldchange (log2FC ~2) seen in PKD1L2, while the lowest was in UBD (log2FC ~1).Despite the absence of time-specific clustering, heatmap analysis indicated the five genes were capable of differentiating between ZIKV-infected and healthy samples (Fig 5A).Moreover, PCA was applied using these five genes, with specific sample clustering identified at 72h and 48h pi (Fig 5B).We further investigated levels of gene expression of each of the five DEGs across time (Fig 5C).Similar levels of gene expression were found for all five genes, with progressive increases observed at 24h pi and 48h pi, respective peaks at 72h pi, followed by dramatically lower expression levels after 5 weeks pi.

Discussion
Both human and animal studies have detected ZIKV in neuronal cells, including progenitor and mature stem cells, as well as astrocytes in vitro.The virus has also been detected in neurosensory tissue, such as the retina and optic nerve [34].Research has shown the alterations viruses cause in cellular pathways create optimal environments for their replication.These strategies make energy resources available for viruses, provide substrate for viral particles and increase the survival of infected cells to consequently ensure viral replication [35].Accordingly, enhancing our understanding of the metabolic alterations induced by viruses may yield useful insights for therapeutic applications in viral diseases.The exact mechanisms (biochemical and cellular) by which the neuropathogenesis of ZIKV happens remains poorly understood [36].
Recent in vitro studies have demonstrated that ZIKV avoids the IFN response (pathway associated with IFN-stimulated response) by destroying STAT2 via proteasomal degradation.The present study conducted transcriptomic analysis on neural cells infected with ZIKV at several post-infection timepoints in order to identify not only differentially expressed genes in comparison to uninfected controls, but also across different time intervals.We further investigated pathways associated with infection at each time interval.In addition to molecular/biochemical alterations caused by ZIKV in host cells on a transcriptomic level, these can also be observed at the proteomic level.Another study investigating the proteome involved in ZIKV infection in human mesenchymal stem cells (hMSC) observed that over 100 proteins were altered by viral infection, many known to be related to different neuropathologies [37].A systematic review of several studies investigating ZIKV-infected organoids identified that viral load was increased after 72 hours regardless of ZIKV viral lineage, a finding consistent with our data that indicates a peak in DEG detection at 72 hours post-infection [37,38].
The present study identified five DEGs shared across all seven-time intervals evaluated: PLA2G2F, TMEM71, PKD1L2, UBD, and TNFAIP3.Another study investigating the effect of ZIKV infection on NSC (neural stem cells) in miRNA signatures and gene expression reported alterations in several signaling pathways due to ZIKV [14].Similarly to the present study, these authors observed increased TMEM71 expression (2.5-, 8.7-, and 35.8-fold increases respectively after 24, 48, and 72 hours of infection) and PKD1L2 (2.8-, 9.6-, and 49.9-fold increases at 24, 48, and 72 after infection, respectively).The authors further reported the dysregulation of several pathways involved in nervous system development and cell cycle, and associated these dysregulations with the formation of microcephaly.Finally, a study conducted in non-human primates infected with ZIKV observed an up-regulation of genes related to interferon, which stands in agreement with the present results [39].
In another study that used NPCs (neuronal precursor cells) to investigate the effect of ZIKV infection on cellular pathways, the authors attempted to establish correlations with clinical aspects of disease (microcephaly, epilepsy, etc.), observing that viral infection activates a variety of inflammatory signaling, similar to our observations after 48h in affected pathways [40].Furthermore, in 2017, a group of researchers investigated the presence of ZIKV in the cerebrospinal fluid (CSF) and lymph nodes (LN) in rhesus monkeys.These authors detected upregulation in the mTOR, proinflammatory, and anti-apoptotic signaling pathways, with the simultaneous downregulation of extracellular matrix signaling pathways.They reported the strong induction of IFN-alpha and other antiviral pathways at 2, 4 and 6 days post-infection, as well as the upregulation of inflammasome components, cytokines and immunomodulatory pathways [41].The data shows a high gene variation in the 72h post-infection.Since these results are based on available data, information regarding viral load is not available.This limitation impairs the identification of an association between viral activity and gene expression variation.
In conclusion, the present findings provide evidence that ZIKV infection in human neuronal cells induced the differential expression of several genes that varied across a range of time intervals post-infection.MDP analysis revealed the 72h pi timepoint as significant for heterogenous gene expression, and the DEGs identified were linked to interleukin and extracellular matrix signaling pathways.Moreover, five DEGs observed across all timepoints were found to be capable of distinguishing between ZIKV-infected and healthy hNSCs.Taken together, those genes could serve as a workbench for future works using RNA-seq data.This methodology allows us to estimate both host and viral RNA expression.This information together would provide useful insights into the viral/host relationship and cell damage process of neuropathogenesis.

Fig 1 .
Fig 1. Gene expression data sampling flow chart.Sampling process showing the query terms in the database, the exclusion criteria and the final number of datasets used in the downstream analysis.https://doi.org/10.1371/journal.pone.0290209.g001

Fig 4 .
Fig 4. Variation analysis of ZIKV: Molecular degree of perturbation analysis in the ZIKV infected samples using healthy samples as baseline.(A) Box Plot of MDP scores over time.Medians were tested using Kruskal-Wallis with Turkey's post-hoc test.(B) Bar plot represent the value observed in each sample.https://doi.org/10.1371/journal.pone.0290209.g004

Fig 5 .
Fig 5. Gene expression levels over time.(A) Heatmap generated from gene expression levels of five DEGs identified at all time points.The colored line (above gene names) reflects the grouping of samples according to gene expression level.The sidebar plot depicts significant log2 fold-change values (FDR<0.05)for each gene across the evaluated timepoints compared to healthy controls.The bottom bar plot reflects the MDP score values of each sample.(B) Gene expression levels of 5 shared DEGs across the evaluated timepoints: Squares indicate median expression levels at each respective timepoint, with bars representing standard error.https://doi.org/10.1371/journal.pone.0290209.g005